#*************         NHANES         *************
#*************      Analyze Code      *************
#************* Logistic 回归模型讲解  *************
#*************                        *************
#### 0.准备好环境 ####
library(tidyverse)
library(gtsummary)
library(tidyr) # drop_na 函数，用于快速去掉 NA
library(survey)
library(haven)


#### 1. 常规Logistics模型 ####
##### 1.1 GLM 函数 #####
# Logit模型结果(是/否)所包含的信息要比定量结果少得多。因此，拟合一个合理的 Logit模型需要的数据比一个相似大小的线性模型要多。
#一个 Logit模型通常至少需要100个观测数据才能符合截距，另外每加预测因子（核心变量+协变量）还需要8-15个观测数据
# Logit 模型中Y变量的“是”或者“否”的占比也会影响到所需要的数据量，
# 通常可以以这个方式进行估计：100+10 * (#X)/(较小占比的Y的结果）
# e.g. #X = 5；Y-1-10%； 100+ 10*5/0.1 = 100+500
  
# 生成 weight 变量，100 个男性 & 100 个女性
n = 100
weight.male <- rnorm(n, 75, 10)
weight.female <- rnorm(n, 55, 5)
weight <- c(weight.male, weight.female)

# 生成性别变量，100 个男性 & 100 个女性
sex.male <- factor(sample(c('male'), n, TRUE))
sex.female <- factor(sample(c('female'), n, TRUE))
sex <- c(sex.male, sex.female)


# 如果要随机生成 100 个男性和女性
# sex <- factor(sample(c('male', 'female'), n, TRUE))

# 生成身高变量
height <- weight + rnorm(n * 2, 0, 5) # 身高变量为 weight + 一个随机的、有点大的扰动项

# 构建简单的Logistic回归模型
# 注意：我们使用性别变量作为 Y，看下 weight、height 能否预测性别
# 要先把 Y 变为 0、1 的值, 更直观，不变的话也不会报错，会自动识别
sex.0.1 <- ifelse(sex == 'male', 1, 0)
data <- data.frame(height = height, weight = weight, sex = sex, sex.0.1 = sex.0.1)

m1 <- glm(
  sex.0.1 ~ weight + height, #
  data = data,
  family = binomial(link = "logit")# 这里也可以直接写 family = "bionomial", 因为 logit 链接是二项式家族的默认方法
)
summary(m1) # 系数输出提供了估计系数及其标准误差，再加上 Wald Z 统计量，即估计系数除以其标准误差。将其与标准的 Normal 分布进行比较，以获得在 Pr (> | z |)列中总结的双尾 p 值。


coefficients(m1)
exp(coefficients(m1)) # weight 前的系数表明，odds 的笔直随着weight每增加1，升高1.42 ，身高越大，Y取1（男性）的概率越大
confint(m1) #与线性回归一样，我们可以得到截距系数和斜率系数的95% 置信区间(为了得到其他水平，改变限制中的水平参数)。
fitted(m1)
AIC(m1) # AIC 越低代表模型的预测效果越好
predict(m1)

##### 1.2 使用 gtsummary 优化输出 #####

m1_tbl_1 <-tbl_regression(m1)
m1_tbl_1

##### 1.3 可视化结果  #####
ggplot(data, aes(x = weight, y = sex.0.1)) +
  geom_jitter(height = 0.05) +
  geom_line(aes(x = weight, y = fitted(m1)), 
            col = "blue") +
  labs(title = "Logistic Regression from Model")


#### 2. 加权 Logistic 模型 ####
##### 2.1 提取数据模块 ##### 
demo.i <- read_xpt("2015-2016/Demographics/demo_i.xpt")#要提前设置好数据存储的路径
colnames(demo.i)

##### 2.2 提取变量  & 去掉 NA 的行 & 衍生变量 ##### 
analyze.sample.data <- demo.i[,c('SEQN',
                                 "RIDAGEYR", "RIAGENDR", "INDFMPIR", "DMDEDUC2",
                                 "WTINT2YR", "SDMVPSU", "SDMVSTRA")]
dim(analyze.sample.data) # 9971    
# View(analyze.sample.data) 

# 一键去掉 NA 的行
analyze.sample.data.drop.na <- drop_na(analyze.sample.data) 
dim(analyze.sample.data.drop.na) # 5082    

# 衍生低收入 vs 非低收入的变量 PIR.factor, 作为后续演示的 Y 变量
analyze.sample.data.drop.na$PIR.factor <- ifelse(analyze.sample.data.drop.na$INDFMPIR < 1.3, 1, 0)


##### 2.3 生成复杂抽样 NHANES_design ##### 
NHANES_design <- svydesign(
  data = analyze.sample.data.drop.na, 
  ids = ~SDMVPSU, 
  strata = ~SDMVSTRA, 
  nest = TRUE, 
  weights = ~WTINT2YR,
  survey.lonely.psu = "adjust") # 可以加上 survey.lonely.psu = "adjust" 避免1个PSU报错

summary(NHANES_design)

##### 2.4 复杂抽样的 Logistics 回归模型 ##### 
# 年龄、教育程度和贫困指数之间的相关性
# family 的选择-binomial，的默认参数是 logit，binomial(link = "logit")

m1 <- svyglm(PIR.factor ~ RIDAGEYR + factor(DMDEDUC2), design = NHANES_design,
             family = binomial)

summary(m1)
# 如果 family 的参数是： binomial，很有可能跳出提示：In eval(family$initialize) : non-integer #successes in a binomial glm!
# 不要紧，可以忽略，glm 只是在指定二项式(和泊松)模型时对原始数据的分布比较谨慎，
# 但不会影响结果，无论如何符合模型。如果你想隐藏这个警告，可以使用 family = quasibinomial 来代替

# quasibinomial 主要的应用场景，若样本观测值变异性过大，即出现了过度离散现象，此时仍使用二项分布假设-binomial就会影响系数检测的显著性。
# 那么补救的方法是使用准二项分布（quasibinomial）。首先要检测样本是否存在过度离散现象，方法是用残差除以残差自由度，若超过1则意味着过度离散。那么将family参数改为quasibinomial。
# quasibinomial(link = "logit") quasibinomial 的默认参数是 logit

m1 <- svyglm(PIR.factor ~ RIDAGEYR + factor(DMDEDUC2), design = NHANES_design,
             family = quasibinomial)
summary(m1)


##### 2.5 使用 gtsummary 优化输出 #####
tbl_regression(m1)

# 显示更多的参数
tbl_regression(m1, exponentiate = TRUE)%>%  # 得到 OR exponentiate model coefficients
  add_global_p() %>%
  add_glance_table(
    include = c(nobs, AIC, BIC))


